clear
est clear


cap mat define migrationfigure=J(105,10,.)
		local row=1

	use "$input/full_geoid_pair_panel.dta"
	

			
	replace migration=migration/(tot_pop/100000)

	//Make sure only include states that have a storm in at least one county at least once in our sample years
	bysort state_fips: egen max=max(storm)
		drop if max!=1
		
	//Create future pooled migration - total migrants in the subsequent five years
	sort geoid_pair move_year
	g post_storm=1 if storm==1|storm[_n-1]==1&geoid_pair==geoid_pair[_n-1]|storm[_n-2]==1&geoid_pair==geoid_pair[_n-2]|storm[_n-3]==1&geoid_pair==geoid_pair[_n-3]|storm[_n-4]==1&geoid_pair==geoid_pair[_n-4]|storm[_n-5]==1&geoid_pair==geoid_pair[_n-5]
	mvencode post_storm, mv(0)
	
	//Collapse to total migration from a county (not looking at paired county to county migration, just all migration out of a county)
	collapse (sum) migration (mean) share_approved totalapprovedihpamount,by(from_geoid move_year storm post_storm state_fips)
	
	//Create IHS transformations
	ihstrans(migration)
	
	
	//Create lags of storms (was there a storm in t-1 years)

			
	label var post_storm "1(0-5 yrs post storm)"
	label var storm "1(Storm year)"
	sort from_geoid move_year
	
	replace totalap=totalap/1000000
	g highdamage=(totalapp>10) if totalapp!=.
	mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
	foreach var in L1 L2 L3 L4 L5{
		label var highdamage_`var' "1(Storm year), `var'"
	}
		
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			
			eststo t1c1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="no lags, net"
			local row=`row'+1
		reghdfe ihs_migration i.highdamage i.highdamage_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L1]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L1]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L2]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L2]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L3]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L3]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L4]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L4]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L5]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L5]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			
			
			lincom 1.highdamage+1.highdamage_L1+1.highdamage_L2+1.highdamage_L3+1.highdamage_L4+1.highdamage_L5
				local beta11=trim("`: display %10.3f r(estimate)'")
				local se11=trim("`: display %10.3f r(se)'")	
				cap mat def migrationfigure[`row',1]=`beta11'
				cap mat def migrationfigure[`row',2]=`se11'
				cap mat def migrationfigure[`row',3]="lags, out"
				local row=`row'+1	
		
		drop highdamage*
		g highdamage=(totalapp>20) if totalapp!=.
		mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
		foreach var in L1 L2 L3 L4 L5{
			label var highdamage_`var' "1(Storm year), `i'"
		}
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)

			eststo t1c2
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="no lags, net"
			local row=`row'+1
		reghdfe ihs_migration i.highdamage i.highdamage_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c2
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L1]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L1]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L2]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L2]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L3]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L3]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L4]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L4]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L5]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L5]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			
			
			lincom 1.highdamage+1.highdamage_L1+1.highdamage_L2+1.highdamage_L3+1.highdamage_L4+1.highdamage_L5
				local beta21=trim("`: display %10.3f r(estimate)'")
				local se21=trim("`: display %10.3f r(se)'")	
				cap mat def migrationfigure[`row',1]=`beta21'
				cap mat def migrationfigure[`row',2]=`se21'
				cap mat def migrationfigure[`row',3]="lags, out"
				local row=`row'+1	
		
		drop highdamage*
		g highdamage=(totalapp>40) if totalapp!=.
		mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
		foreach var in L1 L2 L3 L4 L5{
			label var highdamage_`var' "1(Storm year), `i'"
		}
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)

		eststo t1c3
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="no lags, net"
			local row=`row'+1
		reghdfe ihs_migration i.highdamage i.highdamage_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c3
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L1]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L1]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L2]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L2]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L3]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L3]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L4]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L4]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L5]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L5]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			
			
			lincom 1.highdamage+1.highdamage_L1+1.highdamage_L2+1.highdamage_L3+1.highdamage_L4+1.highdamage_L5
				local beta31=trim("`: display %10.3f r(estimate)'")
				local se31=trim("`: display %10.3f r(se)'")	
				cap mat def migrationfigure[`row',1]=`beta31'
				cap mat def migrationfigure[`row',2]=`se31'
				cap mat def migrationfigure[`row',3]="lags, out"
				local row=`row'+1	
		
		drop highdamage*
		g highdamage=(totalapp>80) if totalapp!=.
		mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
		foreach var in L1 L2 L3 L4 L5{
			label var highdamage_`var' "1(Storm year), `i'"
		}
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)

		eststo t1c4
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="no lags, net"
			local row=`row'+1
		reghdfe ihs_migration i.highdamage i.highdamage_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c4
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L1]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L1]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L2]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L2]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L3]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L3]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L4]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L4]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L5]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L5]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			
			
			lincom 1.highdamage+1.highdamage_L1+1.highdamage_L2+1.highdamage_L3+1.highdamage_L4+1.highdamage_L5
				local beta41=trim("`: display %10.3f r(estimate)'")
				local se41=trim("`: display %10.3f r(se)'")	
				cap mat def migrationfigure[`row',1]=`beta41'
				cap mat def migrationfigure[`row',2]=`se41'
				cap mat def migrationfigure[`row',3]="lags, out"
				local row=`row'+1
		
		
		drop highdamage*
		g highdamage=(totalapp>160) if totalapp!=.
		mvencode highdamage, mv(0) o
		sort from_geoid move_year
		foreach i in 1 2 3 4 5{
				g highdamage_L`i'=highdamage[_n-`i'] if from_geoid==from_geoid[_n-`i'] 
		}	
		foreach var in L1 L2 L3 L4 L5{
			label var highdamage_`var' "1(Storm year), `i'"
		}
		reghdfe ihs_migration i.highdamage, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)

		eststo t1c5
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="no lags, net"
			local row=`row'+1
		reghdfe ihs_migration i.highdamage i.highdamage_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c5
			cap mat def migrationfigure[`row',1]=_b[1.highdamage]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L1]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L1]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L2]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L2]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L3]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L3]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L4]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L4]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',1]=_b[1.highdamage_L5]
			cap mat def migrationfigure[`row',2]=_se[1.highdamage_L5]
			cap mat def migrationfigure[`row',3]="lags, out"
			local row=`row'+1
			
		lincom 1.highdamage+1.highdamage_L1+1.highdamage_L2+1.highdamage_L3+1.highdamage_L4+1.highdamage_L5
				local beta51=trim("`: display %10.3f r(estimate)'")
				local se51=trim("`: display %10.3f r(se)'")	
				cap mat def migrationfigure[`row',1]=`beta51'
				cap mat def migrationfigure[`row',2]=`se51'
				cap mat def migrationfigure[`row',3]="lags, out"
				local row=`row'+1	
				
				
local row=1
local c=1
foreach i in 10 20 40 80 160{
	//Create future return migration panel
	{
			clear
			use "$input/full_geoid_pair_panel.dta"
				
			replace migration=migration/(tot_pop/100000)

			
				
			sort from_geoid to_geoid move_year
			
			keep from_geoid to_geoid move_year migration
			
			rename migration migration_into_from_out_to
			
			rename (from_geoid to_geoid) (to_geoid from_geoid)
			
			save "$input/temp1.dta", replace
			
				clear
			use "$input/full_geoid_pair_panel.dta"
				
			replace migration=migration/(tot_pop/100000)

			
				
			sort from_geoid to_geoid move_year
			
			keep from_geoid to_geoid move_year migration
			
			cap drop _merge
			merge 1:1 from_geoid to_geoid move_year using "$input/temp1.dta"
			mvencode migration_into_from_out_to, mv(0) o
			*drop if _merge!=3
			g paired_net_migration=migration_in-migration
			
			drop _merge
			
			save "$input/paired_net_migration.dta", replace
			clear
	}
	
	use "$input/full_geoid_pair_panel.dta"
		
	replace migration=migration/(tot_pop/100000)

	//Make sure only include states that have a storm in at least one county at least once in our sample years
	bysort state_fips: egen max=max(storm)
		drop if max!=1
	
	replace totalap=totalap/1000000
	g highdamage=(totalapp>`i') if totalapp!=.
	
	keep if highdamage==1 & migration>0 & migration!=.
		
	keep from_geoid to_geoid move_year
	duplicates drop from_geoid to_geoid move_year, force
	g high_damage_receiving_county=1
		save "$input/temp.dta", replace
	clear
	
	use "$input/full_geoid_pair_panel.dta"
		replace migration=migration/(tot_pop/100000)
		bysort state_fips: egen max=max(storm)
			drop if max!=1
		
		replace totalap=totalap/1000000
		g highdamage=(totalapp>`i') if totalapp!=.
		merge 1:1 from_geoid to_geoid move_year using "$input/temp.dta", nogen keep(1 3)
	
		cap drop _merge
		merge 1:1 from_geoid to_geoid move_year using "$input/paired_net_migration.dta"
		mvencode highdamage high_damage_receiving_county, mv(0) o
		drop if _merge==2
		bysort from_geoid to_geoid: egen received_storm_migrations=max(high_damage_receiving_county)
		
		drop if received_storm_migrations!=1
		
			
		//Create IHS transformations
		ihstrans( paired_net_migration)
	
		label var high_damage_receiving_county "1(Storm year)"
		
		sort from_geoid to_geoid move_year
		foreach i in 1 2 3 4 5{
				g high_damage_receiving_county_L`i'=high_damage_receiving_county[_n-`i'] if from_geoid==from_geoid[_n-`i'] & to_geoid==to_geoid[_n-`i']
		}	
		foreach var in L1 L2 L3 L4 L5{
			label var high_damage_receiving_county_`var' "1(Storm year), `i'"
		}
		
		
		
		reghdfe ihs_paired_net_migration i.high_damage_receiving_county, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t1c`c'
			cap mat def migrationfigure[`row',4]=_b[1.high_damage_receiving_county]
			cap mat def migrationfigure[`row',5]=_se[1.high_damage_receiving_county]
			cap mat def migrationfigure[`row',6]="no lags, net"
			local row=`row'+1
		reghdfe ihs_paired_net_migration i.high_damage_receiving_county i.high_damage_receiving_county_L*, absorb(i.from_geoid  i.move_year i.state_fips#c.move_year) vce(cluster state_fips)
			eststo t2c`c'
			cap mat def migrationfigure[`row',4]=_b[1.high_damage_receiving_county]
			cap mat def migrationfigure[`row',5]=_se[1.high_damage_receiving_county]
			cap mat def migrationfigure[`row',6]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',4]=_b[1.high_damage_receiving_county_L1]
			cap mat def migrationfigure[`row',5]=_se[1.high_damage_receiving_county_L1]
			cap mat def migrationfigure[`row',6]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',4]=_b[1.high_damage_receiving_county_L2]
			cap mat def migrationfigure[`row',5]=_se[1.high_damage_receiving_county_L2]
			cap mat def migrationfigure[`row',6]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',4]=_b[1.high_damage_receiving_county_L3]
			cap mat def migrationfigure[`row',5]=_se[1.high_damage_receiving_county_L3]
			cap mat def migrationfigure[`row',6]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',4]=_b[1.high_damage_receiving_county_L4]
			cap mat def migrationfigure[`row',5]=_se[1.high_damage_receiving_county_L4]
			cap mat def migrationfigure[`row',6]="lags, out"
			local row=`row'+1
			cap mat def migrationfigure[`row',4]=_b[1.high_damage_receiving_county_L5]
			cap mat def migrationfigure[`row',5]=_se[1.high_damage_receiving_county_L5]
			cap mat def migrationfigure[`row',6]="lags, out"
			local row=`row'+1	
			
			lincom 1.high_damage_receiving_county+1.high_damage_receiving_county_L1+1.high_damage_receiving_county_L2+1.high_damage_receiving_county_L3+1.high_damage_receiving_county_L4+1.high_damage_receiving_county_L5
				local beta`c'1=trim("`: display %10.3f r(estimate)'")
				local se`c'1=trim("`: display %10.3f r(se)'")	
				cap mat def migrationfigure[`row',4]=`beta`c'1'
				cap mat def migrationfigure[`row',5]=`se`c'1'
				cap mat def migrationfigure[`row',6]="lags, out"
				local row=`row'+1	
			
			local c=`c'+1
	clear	
	
}		
	
	
	
	
	
	
	preserve
		clear
			
		svmat2 migrationfigure

		drop if migrationfigure1==.
		
		keep in 1/8
	
		g n=_n
			replace n=1.8 if n==2
			replace n=2.4 if n==3 
			replace n=3 if n==4
			replace n=3.6 if n==5
			replace n=4.2 if n==6
			replace n=4.8 if n==7
			replace n=5.6 if n==8
		g n2=n+0.25	
		
		g HI95_mean=(1.96*migrationfigure2)+migrationfigure1
		g LI95_mean=migrationfigure1-(1.96*migrationfigure2)

		g HI99_mean=(2.06*migrationfigure2)+migrationfigure1
		g LI99_mean=migrationfigure1-(2.06*migrationfigure2)
		
		g HI95_mean_n=(1.96*migrationfigure5)+migrationfigure4
		g LI95_mean_n=migrationfigure4-(1.96*migrationfigure5)

		g HI99_mean_n=(2.06*migrationfigure5)+migrationfigure4
		g LI99_mean_n=migrationfigure4-(2.06*migrationfigure5)
		
		set scheme plotplain
		
		rename migrationfigure1 mean_beta
		rename migrationfigure4 mean_beta_n
		
		replace mean_beta=mean_beta*100 
		replace HI99_mean=HI99_mean*100 
		replace LI99_mean=LI99_mean*100
		
		replace mean_beta_n=mean_beta_n*100 
		replace HI99_mean_n=HI99_mean_n*100 
		replace LI99_mean_n=LI99_mean_n*100
		
		replace n=n-1
		replace n2=n2-1
	
		twoway 	(bar mean_beta_n n2  if n<8, barwidth(0.25) color("252 89 16%50")) ///
					(rspike HI99_mean_n LI99_mean_n n2 if n<7, lc(gs8))  ///
					(bar mean_beta n  if n<8, barwidth(0.25) color("0 45 114%50")) ///
					(rspike HI99_mean LI99_mean n if n<8, lc(gs8)  ///
					yti("Percentage change in migration",size(vlarge)) ylabel(, nogrid labsize(large)) xlabel(0 "Storm year" 0.8 "Storm year" 1.4 "Storm year, L1" 2 "Storm year, L2" 2.6 "Storm year, L3" 3.2 "Storm year, L4" 3.8 "Storm year, L5" 4.6 "Sum of lags",nogrid labsize(large)) ysize(1.5)  xsize(5) xtitle("") legend(off) yline(0) xline(0.4) xline(4.2))
					
					graph export "$figures/fig_2a.png", as(png) replace
					
	restore	
	
	
	preserve
		clear
			
		svmat2 migrationfigure

		drop if migrationfigure1==.
		
		keep in 9/16
	
		g n=_n
			replace n=1.8 if n==2
			replace n=2.4 if n==3 
			replace n=3 if n==4
			replace n=3.6 if n==5
			replace n=4.2 if n==6
			replace n=4.8 if n==7
			replace n=5.6 if n==8
		g n2=n+0.25	
		
			g HI95_mean=(1.96*migrationfigure2)+migrationfigure1
		g LI95_mean=migrationfigure1-(1.96*migrationfigure2)

		g HI99_mean=(2.06*migrationfigure2)+migrationfigure1
		g LI99_mean=migrationfigure1-(2.06*migrationfigure2)
		
		g HI95_mean_n=(1.96*migrationfigure5)+migrationfigure4
		g LI95_mean_n=migrationfigure4-(1.96*migrationfigure5)

		g HI99_mean_n=(2.06*migrationfigure5)+migrationfigure4
		g LI99_mean_n=migrationfigure4-(2.06*migrationfigure5)
		
		set scheme plotplain
		
		rename migrationfigure1 mean_beta
		rename migrationfigure4 mean_beta_n
		
		replace mean_beta=mean_beta*100 
		replace HI99_mean=HI99_mean*100 
		replace LI99_mean=LI99_mean*100
		
		replace mean_beta_n=mean_beta_n*100 
		replace HI99_mean_n=HI99_mean_n*100 
		replace LI99_mean_n=LI99_mean_n*100
		
		replace n=n-1
		replace n2=n2-1
	
		twoway 	(bar mean_beta_n n2  if n<8, barwidth(0.25) color("252 89 16%50")) ///
					(rspike HI99_mean_n LI99_mean_n n2 if n<7, lc(gs8))  ///
					(bar mean_beta n  if n<8, barwidth(0.25) color("0 45 114%50")) ///
					(rspike HI99_mean LI99_mean n if n<8, lc(gs8)  ///
					yti("   ",size(vlarge)) ylabel(, nogrid labsize(large)) xlabel(0 "Storm year" 0.8 "Storm year" 1.4 "Storm year, L1" 2 "Storm year, L2" 2.6 "Storm year, L3" 3.2 "Storm year, L4" 3.8 "Storm year, L5" 4.6 "Sum of lags",nogrid labsize(large)) ysize(1.5)  xsize(5) xtitle("") legend(off) yline(0) xline(0.4) xline(4.2))
					
					graph export "$figures/fig_2b.png", as(png) replace
					
					
	restore	
	
	preserve
		clear
			
		svmat2 migrationfigure

		drop if migrationfigure1==.
		
		keep in 17/24
	
		g n=_n
			replace n=1.8 if n==2
			replace n=2.4 if n==3 
			replace n=3 if n==4
			replace n=3.6 if n==5
			replace n=4.2 if n==6
			replace n=4.8 if n==7
			replace n=5.6 if n==8
		g n2=n+0.25	
		
			g HI95_mean=(1.96*migrationfigure2)+migrationfigure1
		g LI95_mean=migrationfigure1-(1.96*migrationfigure2)

		g HI99_mean=(2.06*migrationfigure2)+migrationfigure1
		g LI99_mean=migrationfigure1-(2.06*migrationfigure2)
		
		g HI95_mean_n=(1.96*migrationfigure5)+migrationfigure4
		g LI95_mean_n=migrationfigure4-(1.96*migrationfigure5)

		g HI99_mean_n=(2.06*migrationfigure5)+migrationfigure4
		g LI99_mean_n=migrationfigure4-(2.06*migrationfigure5)
		
		set scheme plotplain
		
		rename migrationfigure1 mean_beta
		rename migrationfigure4 mean_beta_n
		
		replace mean_beta=mean_beta*100 
		replace HI99_mean=HI99_mean*100 
		replace LI99_mean=LI99_mean*100
		
		replace mean_beta_n=mean_beta_n*100 
		replace HI99_mean_n=HI99_mean_n*100 
		replace LI99_mean_n=LI99_mean_n*100
		
		replace n=n-1
		replace n2=n2-1
	
		twoway 	(bar mean_beta_n n2  if n<8, barwidth(0.25) color("252 89 16%50")) ///
					(rspike HI99_mean_n LI99_mean_n n2 if n<7, lc(gs8))  ///
					(bar mean_beta n  if n<8, barwidth(0.25) color("0 45 114%50")) ///
					(rspike HI99_mean LI99_mean n if n<8, lc(gs8)  ///
					yti("   ",size(vlarge)) ylabel(, nogrid labsize(large)) xlabel(0 "Storm year" 0.8 "Storm year" 1.4 "Storm year, L1" 2 "Storm year, L2" 2.6 "Storm year, L3" 3.2 "Storm year, L4" 3.8 "Storm year, L5" 4.6 "Sum of lags",nogrid labsize(large)) ysize(1.5)  xsize(5) xtitle("") legend(off) yline(0) xline(0.4) xline(4.2))
					
					graph export "$figures/fig_2c.png", as(png) replace
					
					
	restore	
	
	preserve
		clear
			
		svmat2 migrationfigure

		drop if migrationfigure1==.
		
		keep in 25/32
	
		g n=_n
			replace n=1.8 if n==2
			replace n=2.4 if n==3 
			replace n=3 if n==4
			replace n=3.6 if n==5
			replace n=4.2 if n==6
			replace n=4.8 if n==7
			replace n=5.6 if n==8
		g n2=n+0.25	
		
			g HI95_mean=(1.96*migrationfigure2)+migrationfigure1
		g LI95_mean=migrationfigure1-(1.96*migrationfigure2)

		g HI99_mean=(2.06*migrationfigure2)+migrationfigure1
		g LI99_mean=migrationfigure1-(2.06*migrationfigure2)
		
		g HI95_mean_n=(1.96*migrationfigure5)+migrationfigure4
		g LI95_mean_n=migrationfigure4-(1.96*migrationfigure5)

		g HI99_mean_n=(2.06*migrationfigure5)+migrationfigure4
		g LI99_mean_n=migrationfigure4-(2.06*migrationfigure5)
		
		set scheme plotplain
		
		rename migrationfigure1 mean_beta
		rename migrationfigure4 mean_beta_n
		
		replace mean_beta=mean_beta*100 
		replace HI99_mean=HI99_mean*100 
		replace LI99_mean=LI99_mean*100
		
		replace mean_beta_n=mean_beta_n*100 
		replace HI99_mean_n=HI99_mean_n*100 
		replace LI99_mean_n=LI99_mean_n*100
		
		replace n=n-1
		replace n2=n2-1
	
		twoway 	(bar mean_beta_n n2  if n<8, barwidth(0.25) color("252 89 16%50")) ///
					(rspike HI99_mean_n LI99_mean_n n2 if n<7, lc(gs8))  ///
					(bar mean_beta n  if n<8, barwidth(0.25) color("0 45 114%50")) ///
					(rspike HI99_mean LI99_mean n if n<8, lc(gs8)  ///
					yti("   ",size(vlarge)) ylabel(, nogrid labsize(large)) xlabel(0 "Storm year" 0.8 "Storm year" 1.4 "Storm year, L1" 2 "Storm year, L2" 2.6 "Storm year, L3" 3.2 "Storm year, L4" 3.8 "Storm year, L5" 4.6 "Sum of lags",nogrid labsize(large)) ysize(1.5)  xsize(5) xtitle("") legend(off) yline(0) xline(0.4) xline(4.2))
					
					graph export "$figures/fig_2d.png", as(png) replace
					
					
	restore	
	
	preserve
		clear
			
		svmat2 migrationfigure

		drop if migrationfigure1==.
		
		keep in 33/40
	
		g n=_n
			replace n=1.8 if n==2
			replace n=2.4 if n==3 
			replace n=3 if n==4
			replace n=3.6 if n==5
			replace n=4.2 if n==6
			replace n=4.8 if n==7
			replace n=5.6 if n==8
		g n2=n+0.25	
	
			g HI95_mean=(1.96*migrationfigure2)+migrationfigure1
		g LI95_mean=migrationfigure1-(1.96*migrationfigure2)

		g HI99_mean=(2.06*migrationfigure2)+migrationfigure1
		g LI99_mean=migrationfigure1-(2.06*migrationfigure2)
		
		g HI95_mean_n=(1.96*migrationfigure5)+migrationfigure4
		g LI95_mean_n=migrationfigure4-(1.96*migrationfigure5)

		g HI99_mean_n=(2.06*migrationfigure5)+migrationfigure4
		g LI99_mean_n=migrationfigure4-(2.06*migrationfigure5)
		
		set scheme plotplain
		
		rename migrationfigure1 mean_beta
		rename migrationfigure4 mean_beta_n
		
		replace mean_beta=mean_beta*100 
		replace HI99_mean=HI99_mean*100 
		replace LI99_mean=LI99_mean*100
		
		replace mean_beta_n=mean_beta_n*100 
		replace HI99_mean_n=HI99_mean_n*100 
		replace LI99_mean_n=LI99_mean_n*100
		
		replace n=n-1
		replace n2=n2-1
	
		twoway 	(bar mean_beta_n n2  if n<8, barwidth(0.25) color("252 89 16%50")) ///
					(rspike HI99_mean_n LI99_mean_n n2 if n<7, lc(gs8))  ///
					(bar mean_beta n  if n<8, barwidth(0.25) color("0 45 114%50")) ///
					(rspike HI99_mean LI99_mean n if n<8, lc(gs8)  ///
					yti("   ",size(vlarge)) ylabel(, nogrid labsize(large)) xlabel(0 "Storm year" 0.8 "Storm year" 1.4 "Storm year, L1" 2 "Storm year, L2" 2.6 "Storm year, L3" 3.2 "Storm year, L4" 3.8 "Storm year, L5" 4.6 "Sum of lags",nogrid labsize(large)) ysize(1.5)  xsize(5) xtitle("") legend(pos(6) col(2) order(1 3) label(1 "Net-migration") label(2 "Out-migration")) yline(0) xline(0.4) xline(4.2))
					
					graph export "$figures/fig_2e.png", as(png) replace
					
					
	restore	
			
			
			